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Abstract 


In this paper we are concerned with high- accuracy quadrature method solutions of nonlinear 
Fredholm integral equations of the form y(x) = r(z)+/o $(*, !/{*))&, 0 < x < 1, where the 

kernel function g(x,t) is continuous, but its partial derivatives have finite jump discontinuities 
across x — t. Such integral equations arise, e.g., when one applies Green’s function techniques 
to nonlinear two-point boundary value problems of the form y"(x) = f(x,y(x)) y 0 < x < 1, 
with y(0) = yo and y(l) = yi, or other linear boundary conditions. A quadrature method 
that is especially suitable and that has been employed for such equations is one based on 
the trapezoidal rule that has a low accuracy. By analyzing the corresponding Euler-Maclaurin 
expansion, we derive suitable correction terms that we add to the trapezoidal rule, thus obtaining 
new numerical quadrature formulas of arbitrarily high accuracy that we also use in defining 
quadrature methods for the integral equations above. We prove an existence and uniqueness 
theorem for the quadrature method solutions, and show that their accuracy is the same as that 
of the underlying quadrature formula. The solution of the nonlinear systems resulting from the 
quadrature methods is achieved through successive approximations whose convergence is also 
proved. The results are demonstrated with numerical examples. 
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1 Introduction 


Consider the Fredholm integral equation of the second kind of the form 

y(x) = r(x) + [ g(x,t)F(t,y(t))dt, 0 < x < 1, (1.1) 

JO 

where the function F(t , w ) is assumed to be nonlinear in w, in general. Let M be a nonnegative 
integer and assume the following: 

(i) r € C M (I), where I = [0, 1]. 

(ii) g e C(ft), where Q = I X I. If M > 1, then the partial derivatives g(x,t) = g j:k {x ,t) 

with j > 0, k > 0, and 1 < j + k < M , are all in PC(fl). By this we mean that they 
are continuous in each of the two halves S_ = {(x,t) : 0 < t < x < 1} and S+ = {(x,t) : 
0 < x < t < 1} of Q, but they are discontinuous across the diagonal S + fi S_ of Q, i.e., 
across x = t, where they have finite jump discontinuities. For future reference let us define 
S k (x) = y 0ife (x,x+ ) — 9o,k{ x i x ~)i k = 1, 2, M. By the assumptions above, S k {x) are 
continuous on I and thus bounded there. 

(iii) F{t, w) £ C{ A) and also F 0>1 {t, w) = £F{t, w) € C(A), where A = / x J with J — [Pi, R 2 ] 
for some R x and R 2 that can be finite or infinite. For M > 3 we also assume that F j>k {t , w) = 

pit, w ), with j + k < M -2, are all in C(A). (Starting with our discussion of improved 
quadrature methods in Section 3, we will assume this with j + k < M for M > 1.) 

Thus, for each value of M , the assumptions in (i)-(iii) contain those for lower values. In 
particular, we have r € C{I), g € C(fi), and F,F 0 ,i € C(A), for any M > 0. These minimal 
smoothness conditions on r, y, and F, along with other conditions not pertaining to smoothness, 
are sufficient to guarantee the existence and uniqueness of (i) a continuous solution y(x) of (1.1), 
cf. Theorem 2.1, and (ii) a quadrature method (approximate) solution of (1.1), cf. Theorem 5.1. 
Theorem 2.1 in the next section states, furthermore, that y(x) € C M (I) for each M > 0 under the 
conditions of (i)— (iii). In particular, y € C°°(7) when M = oo. 

Integral equations of the kind described in this introduction arise, for example, when one applies 
Green’s function techniques to two-point boundary value problems (BVP’s) governed by nonlinear 
ordinary differential equations (ODE’s). See, e.g., Courant and Hilbert [CH], Morse and Feshbach 
[MF], Keller [K], and Pennline [PI], [P2], and [P3]. 
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To illustrate this point let us consider 


y"(x) = f(x,y(x)), o < * < 1 , 


with the inhomogeneous boundary conditions (BC’s), 


( 1 . 2 ) 


y(0) = y o and y(l) = y x . 


(1.3) 


As is shown in [K], (1.2)— (1.3) can be converted into the Fredholm integral equation of the second 
kind 


y{x) = r(x) + J g(x, t)[k 2 y(t ) - f(t, y{t))]dt, 0 < x < 1, 


(1.4) 


where 


g{x,t) 


l 

A: sinh k 


sinh fcx sinh fc(l — t), 0 < x < t 
sinh k( 1 — x) sinh kt, t < x < 1 


(1.5) 


and 


/ v _ y 0 sinh fc( 1 - x) + yi sinh kx 

sinh k ' ' 

Here k > 0 is a free parameter chosen to guarantee the convergence to the solution y(x) of the 
sequence of successive approximations {y( m *(x)}{£_ 0 obtained as in 


y( m +i) (3.) = r(x) + £ g(x,t)[k 2 yM(t) - f(t,yM(t)))dt , m = 0, 1, ... , (1.7) 


with chosen suitably. 

A standard procedure for solving (1.1) numerically is the quadrature method; see, e.g., Baker 
[B, p. 686]. In this method we start with a numerical quadrature formula Ijv[4>] = 
for the integral /J <j>(t)dt. Here 0 < x 0 < x x < • • • < x N < 1. Next, we replace the integral 
Jo g(x , t)F(t , y(t))dt by the corresponding *)F]. Finally, we collocate the resulting equation 

at the abscissas z, , i — 0, 1, N , to obtain the nonlinear system of equations 

N 

yi = r(x.) + J2 a j9i x i , Xj)F(xj,yj), i = 0, 1, ..., N, (1.8) 

i= o 

where, for each i, j/j is the approximation to 

Subsequently, this system may be solved, e.g., by successive approximations as in 


yl 0) = y {0) (xi), i = o,i, 


N 


y< m+1) = r(x,) + a i9{Xi, x i )F(x j , yj m) ), i = 0, 1, ..., JV; to = 0,1,.... (1.9) 

j=0 
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One can also use Newton’s method for solving the system in (1.8), but this requires the computation 
of the Jacobian matrix and the solution of a linear system of N + 1 equations at each iteration, 
which may make the solution very expensive computationally. See, e.g., [K] and [B]. We shall come 
back to this subject in Section 8, where we will discuss other options as well. 

In general, the accuracy of the y, in (1.8) is that provided by the numerical quadrature formula 
I N [g(x, -)F]> subject to the condition that g(x,t)F{t, y(t)) is sufficiently smooth for t £ I. For the 
case considered in this work, however, g(x,t)F(t, y(t)) is not continuously differentiable for t £ /, 
but only continuous there. This is so since g 01 (x,t) = jfcg(x,t) has a (finite) jump discontinuity 
for t = x. Therefore, we cannot expect to obtain a high-accuracy numerical solution by using 
a high-accuracy numerical quadrature formula such as a Gaussian formula. For this reason, the 
trapezoidal rule that has a low accuracy of 0(N~ 2 ) has been used in previous work, see [K]. 

When the approach above, with In taken as the trapezoidal rule, is applied to the integral 
equation (1.4)— (1.6) , the resulting y,- have errors of order 0{N ~ 2 ) as shown in [K], provided that 
y £ C 2 (I) and {y (m) (x)}~ =0 defined by (1.7) is a contractive sequence. The same approach was 
used also in [P1]-[P3]. 

In the present work we propose to improve the accuracy from 0(N~ 2 ) to 0(N~ 2p ) for arbitrary 
integers p > 2, by replacing the trapezoidal rule by “numerical quadrature formulas” that have 
higher accuracy in the presence of the nonsmooth kernels g(x,t) that we consider here. Specifically, 
these formulas are obtained by adding suitable correction terms to the trapezoidal rule approxima- 
tions at the endpoints t = 0 and t = 1 and also at t = x, the point where g{x,t) fails to be smooth. 
These terms are derived from a careful analysis of the Euler-Maclaurin expansion associated with 
the error in the trapezoidal rule. Due to the nature of the correction terms, what we obtain are 
not real numerical quadrature formulas in the sense described in the paragraph following (1.7). 

An important point that will be seen later is that given N, the amount of computational work 
per iteration is practically independent of the order of accuracy N~ 2p of the quadrature formula 
used. This means we can increase the order of accuracy by keeping the cost per iteration almost 
the same. 

An approach similar in spirit to the one here was taken by Sidi and Israeli [SI] in the quadrature 
method solution of periodic Fredholm integral equations with weakly singular kernels that have 
algebraic/logarithmic singularities along the line t — x. In [SI] too the Euler-Maclaurin expansion 
of the trapezoidal rule plays a crucial role in the development of new numerical quadrature formulas 
of high-order accuracy. Only there the periodic nature of the kernel and the solution enables one 
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to propose extrapolated (Romberg-type) formulas to replace the trapezoidal rule. In the present 
case, however, we do not have any periodicity either in the kernel or in the solution, and, therefore, 
we cannot use extrapolated integration formulas. Instead, we use corrected formulas to replace the 
trapezoidal rule. Thus, the approach, methods, and results of the present work are quite different 
from those of [SI]. 

The existence and uniqueness of the solution to the nonlinear system in (1.8) as well as the 
solution to the integral equation in (1.1) has been discussed in [K, Chap. 4] in the context of 
two-point BVP’s described above. Keller’s results are obtained under the condition that F 0) i(£,w) 
is continuous and bounded for tel and for all w E (-oo,+oo). This is a very severe restriction 
on F, however. Most problems of engineering interest do not satisfy this restriction. In many 
applications physical considerations lead one to conclude that the solution is restricted to some 
finite interval. This suggests that it may be feasible to state existence and uniqueness theorems in 
which F 0> i(ri w) is continuous and, therefore, also bounded for t E / and w E J = [Fi, F 2 ] for some 
finite R x and F 2 , the solution satisfying y(x) E J for x E I as well. This view is taken in the series 
of papers by Pennline, who establishes several existence and uniqueness theorems in the context of 
two-point BVP’s. Pennline also shows how these theorems apply to various problems that arise in 
certain engineering applications. R 1 and F 2 are assumed finite also in the present paper. 

Both in [K] and in [P1]-[P3] the existence and uniqueness of the solution to (1.1) is proved by 
establishing that a sequence {y< m )(x)}^_ 0 °f successive approximations from (1.7) contracts and 
thus converges to the solution y(x) of (1.1) uniformly on I. Before this can be done, however, 
one has to show that if the initial approximation y^(x) satisfies y^°^(x) E J for x E /, where 
J is the finite interval mentioned in the previous paragraph, then so do all the other y( m \x). 
(This is not necessary when J is (— 00 , 00 ).) When analyzing the existence and uniqueness of 
the numerical solution defined by the quadrature methods in (1.8) one would like to adopt the 
same approach. That is to say, we would like to be able to show first that the successive ap- 
proximations y^ in (1.9) satisfy y- m) E J, i = 0,1,. ..,7V, for all m = 1,2,..., and use this 
to establish that {y- m \ i = 0, 1, ..., N}™ =0 contracts and thus has a limit {y,, i — 0,1 ,...,V} 
that is the unique solution to (1.8). Although y (0 )(z) E J may imply that y^(x) E J, y,- 0 ^ = 
y w {x { ) € J may not guarantee that E J, due to the error in the numerical quadrature formula 

f° r /o 1 5(^*1 Similarly, the yj 2) that are obtained 

from the y- 1 ^ and the subsequent y- m) may not all lie in J. In short, the analysis of the nonlinear 
system in (1.8) seems to become rather complicated when J is a finite interval. Simply, the con- 
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ditions r € C(I), g € C(Si), and F, F 0 , i € C( A), for which we are able to state an existence and 
uniqueness theorem for the solution of (1.1) do not seem to suffice for a corresponding theorem for 
the approximate solution defined by (1.8). 

In this paper we consider this problem in detail and prove an existence and uniqueness theorem 
for the numerical solution by extending the condition F,Fo,i € C*(A) slightly to read F,Fo t i € 
C(A'), where A' = / x J', where J' = [Ri~t), R 2 + rj] D J, i? > 0 being arbitrarily small. A useful 
feature of our proof technique is the use of the modulus of continuity in many places. This enables 
us to carry out the analysis without resorting to the t-5 formalism that would have to be used 
otherwise. We believe that the idea of employing the modulus of continuity may be applicable in 
other problems of numerical analysis as well. 

For an existence and uniqueness theorem under assumptions that are of a different nature, see 
[B, pp. 689-691]. 

. The plan of this paper is as follows: 

In the next section we consider the convergence of the method of successive approximations for 
(1.1), and state an existence and uniqueness theorem for the solution of (1.1) that relies on successive 
approximations. We also derive an equicontinuity result for the successive approximations y {m) {x) 
that we use later. 

In Section 3 we derive our higher-order “numerical quadrature formulas” that we use in the 
quadrature method by correcting the trapezoidal rule appropriately. In Section 4 we derive error 
bounds for the trapezoidal rule and its modifications that are expressed in terms of moduli of 
continuity and thus are uniform in the a:,-. These bounds form an essential part of the analyses 
given in Sections 5 and 6. With the slight extension F, F 0 ,i € C'(A') that we mentioned above, 
in Section 5 we prove the convergence of the sequences {yj ^}m=o f° r thereby establishing 

the existence and uniqueness of the numerical solution y* , i = 0, as well. With the same 
extension, in Section 6 we analyze the errors in the y f as functions of N. We do this analysis both 
for the trapezoidal rule and for its modifications. We give uniform bounds on |y f - y(x<)| for all 
M > 0. In particular, the bounds for M = 0 and M = 1 are of forms not encountered before. One 
of the conclusions that can be drawn from this analysis is that if y € C (7), AT ^ 1, then by using 
the appropriate modified trapezoidal rule we can achieve an error of order h M , where h = 1/N. 
Finally, in Section 7 we illustrate the new quadrature method and the accompanying theory with 
specific nonlinear two-point BVP’s. 
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As far as is known to us, the quadrature methods proposed in this work and their accompanying 
theory on existence, uniqueness, and convergence of numerical solutions have not been published 
elsewhere previously. 

We mentioned above that nonlinear two-point BVP’s can be formulated as Fredholm integral 
equations of the second kind of the type treated in this work. Thus, the methods of this work 
can also be used for solving numerically two-point BVP’s. Here it can be argued that solving the 
associated ODE’s by finite differences may be less expensive than solving the corresponding integral 
equations as the difference equations that are formed have banded Jacobian matrices and hence may 
be solved efficiently by Newton’s method. The band size increases as the accuracy of the solution 
is increased, however. It is also known that the finite difference approach has great difficulties in 
treating BVP’s with solutions y(x) that vary rapidly on (0, 1), which may occur, for example, in the 
form of very thin boundary layers. The integral equation approach does not seem to have problems 
in producing numerical solutions in a stable manner also for such BVP’s. In addition, the accuracy 
of the numerical solution of the integral equation approach can be increased arbitrarily, practically 
with no extra computational cost. See, e.g., Example 2 in Section 7. Finally, in case y'(0) or y'(l) 
or both are present in the boundary conditions, they have to be discretized with suitable accuracy 
when solving the ODE’s, whereas in the integral equation approach boundary conditions are built 
right into the associated integral equations and require no discretization. 

2 Existence and Uniqueness for Solution of (1.1) 

Let us pick a function j/°)(a:) and generate the functions m = 1,2,..., by the method of 

successive approximations as in 

y {m+1) (x) = r(x)+ [ g(x,t)F(t,y {m) (t))dt, m- 0,1,.... (2.1) 

Jo 

The following theorem gives a set of sufficient conditions for {j/^ rn H a: )}m=o converge, estab- 
lishing the existence and uniqueness of a continuous solution to (1.1) at the same time. 

Theorem 2.1 Denote 

V[u](x) = r(z)+ [ g{x,t)F{t,u{t))dt , (2.2) 

Jo 

and assume that 

u(x) £ J for x € I implies \P[u](:r) € J for x £ I. (2.3) 
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Assume also that r e C(I), g € C(fi), and F, F 0 ,i 6 C( A). Denote the operator L^-norm of g{x, t ) 
on Q and the L^-norm ofF 0A (t,w) on A by [yj and ||F 0)1 ||, respectively, i.e., 


M = 

r 1 

max / \g(x,t)\dt and ||F 0i i|| = max |F 0 ,i(t, ru)|. 
x€l Jo (:,«;)€ A 

(2.4) 

Then , provided that 

f* = M ll^o, ill <i 

(2.5) 

and 

y (0) e C(I) and y {0] {x) e J for x£ I, 

(2.6) 

the following hold: 




(i) yW e C{I) and y (x) G J for x € /, m = 1,2, ... . 

(ii) {y (m) ( x )}m=o converges uniformly on I to a function y(x) such that y € C(I) and y(x) € J 
for 16 /. 

(Hi) y(x) is the unique solution of (1-1)- 

(iv) If, in addition, r{x), g{x,t), and F(i, te) are as described in (i)-(iii) of the first paragraph of 
Section 1 with arbitrary M, then y € C M (I). 

The proof of parts (i)— (iii) of this theorem are almost identical to that of Theorem 4.1.2 in [K, 
pp. 108-109], provided suitable additions and modifications are made in the latter. 

The result of part (iv) can be verified by splitting the integral in (1.1) into the sum / Q r + /j , 
and then differentiating under the integral sign and using induction on M. (The case M — 0 is 
already covered in parts (i)-(iii).) In the course of the proof it also becomes clear that only those 
g jik (x,t) for which j > 1 and j + k < M - 1 and g M fi{x,t) are required to be in PC{Q ) for M > 1 . 

Our next theorem essentially states that, under the conditions of Theorem 2.1, the sequence 
{t/ m )(x )}~ =0 is equicontinuous on /. We state it in terms of the moduli of continuity of r and y (m) 
on I and of g on ft. For the sake of completeness we give the precise definition of this concept. 

Definition: Let Y$ = F(^,...,^„) be defined on a subset X of R n . Then its modulus of 
continuity u>y on X is defined as 

u Y {h) =sup{|Y(|j - Y(£)i € X, |& - £| < h, i = (2.7) 
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It is known that if X is a compact set and Y(£) is continuous on X , therefore, uniformly 
continuous there, wy(/i) — ► 0 monotonically as h — ► 0. We refer the reader to Cheney [C] for this 
and other details on moduli of continuity. 

Theorem 2.2 Define 

11*11 = max |F(t,u;)|, (2.8) 

(t,tu)€ A 

and let cj r , u g , and u y { m ) denote the moduli of continuity of r(x) on I , of g(x,t) on Q f and of 
t/ m )(x) on I , respectively. Then , for any M > 0, we have 

u? y (m)(h) < u* r (h) + \\F\\u>g(h), m = 1,2, (2.9) 

and thus uj y (n)(h) — > 0 as h — > 0 uniformly in m. 

Proof. From (2.1) we have for m = 1,2, ..., 

y {m) (x ) - y <m) (x') = [r(x) - r(x')] + f [y(x, t) - g{x', t)]F(t, y (m_1) (t))dt. (2.10) 

«/o 

The result in (2.9) now follows by taking absolute values on both sides of (2.10) and invoking (2.7) 
and (2.8) along with the result that y^(x) 6 J for x € I. The rest follows from the fact that 
r 6 C(J), y {m) € C(/), m = 0, 1, ..., and g € C(Q). □ 

3 Derivation of the Improved Quadrature Formulas 

Let us denote <f>(t) = g(x, t)F(t, y(t)) with x being held fixed . Let us also assume that, in case 
M > 1, Fj fk (t, w) € C( A) for j + A; < M (instead of j + k < M - 2 for M > 3). Here y(x) is the 
unique solution of (1.1) in the sense of Theorem 2.1, i.e., y E C M (I) and y(x) E J when x E L 
Thus, we are assuming that the conditions (i)-(iii) of Section 1 and the conditions (2.3) and (2.5) 
of Theorem 2.1 are satisfied. We will retain all these assumptions throughout the remainder of this 
work. We conclude that <j>(t) is continuous for t E /, but not continuously differentiable. We also 
conclude that <p^{t), k = 1, M, are continuous in each of the intervals [0, x\ and [x, 1], but have 
finite jump discontinuities at t = x when x E (0, 1). 

Let h = 1/N, where N is some positive integer, and let x t * = i/i, i = 0, 1, ..., N. Assume now 
that the point x mentioned in the previous paragraph is equal to x, for some fixed i. Let us consider 
the trapezoidal rule approximations T-(h) for f* <f>(t)dt and T+(h ) for J 1 <f>(t)dt that are given by 

F.(/ l ) = /i^Xx i ) and T + (h) = hJ2"<f>(x j ), (3.1) 

J=0 j=i 


10 



where = ^(/x r +/x,)+ £ fij if r < a, and = 0 if r = s. Obviously, 

j=r 2 i=r+l j=r 

T_(h) + T + (h) = *£>(*;) = T(h), (3.2) 

j= o 

where T(h ) is the trapezoidal rule approximation for 4>{t)dt. 

3.1 Euler-Maclaurin Expansions for T(h ) 

Let us first consider the cases i = 1, iV — 1. For each such case x = x,- € (0, 1), and we have the 
following (Euler-Maclaurin) expansions for T_(/i) and T+(h ): 

/ r <f>{t)dt = T.{h) - £ TTHT k (25_1> (^-) " ^ (2,_1) (0)] ^ + E { p -\h- x), 

j o 1= i 

££ _) (/i;x) = -^^ 2P) (^)^ P - for some £- 6 (°’ x )’ ( 3 - 3 ) 

and 

t' m* = 2. (ft) -E-ih U (2 — , (l) - ^-‘»(*+)] !>” + E ( »(h; x), 

Jx 5=1 \ ZS )' 

E( +) (h; x) = -(1 - x)^^4> {2p) (Z+)h 2p , for some f+ € (*, 1). (3.4) 

In (3.3) and (3.4) B r are the Bernoulli numbers. Combining (3.3) and (3.4), we have 
fj(t)dt=T(h ) - 

+ £ 7 ^ U ! -‘>(X+) - *-*>(.-)] h-“ + E r (h-,x), 

E,(ft;x) = B<-)(fc; x) + £<«(/>; x) = [x^’K-) + (1 - x)<»< 2rt (« + )] (3.5) 

Let us now turn to the cases i = 0 and i = N. When x = x 0 = 0, we have simply T(h ) = T + {h). 
Therefore, the Euler-Maclaurin expansion of T{h ) now is given by (3.4) with x there replaced by 
0. Similarly, when x = x N = 1, T(/i) = T_(/i), and the Euler-Maclaurin expansion of T{h) is given 
by (3.3) with x there replaced by 1. In other words, T(/i) satisfies (3.5) also when x = Xo = 0 and 
x = x w = 1, with the second summation involving [<?!> (23-1) (x-|-) - </> (25-1) (x-)] being absent from 
(3.5) in both cases. 

Notice that E p (h\ x) = 0{h 2p ) as h ->■ 0 uniformly in x = x { , i - 0, 1, ..., N , and in N . Also in 
(3.3)-(3.5) we have assumed that M > 2p (p > 1). We shall make this assumption throughout the 
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remainder of this section even though it does not cover all possible cases. (We will consider the 
remaining cases following Theorem 6.1 in Section 6.) 

The combined Euler-Maclaurin expansion of (3.5) guides the derivation of the improved numer- 
ical quadrature formulas below. For a discussion of the Euler-Maclaurin expansion see, e.g., Davis 
and Rabinowitz [DR]. 

3.2 Corrections to the Trapezoidal Rule 

It is clear from the Euler-Maclaurin expansions given in (3.3)-(3.5) that <f>(t)dt — T(h) = 0(h 2 ) 
as h -» 0, for all x = x,, i = 0, 1, While for x = x 0 = 0 and x = x N = 1 this result is 

immediate, for x = x, 6 (0, 1) it comes somewhat as a surprise, as <f>'(t) is not continuous on [0, 1] 
for such values of x. We now aim at improving the accuracy of T(h) by taking its Euler-Maclaurin 
expansion into account. 

To motivate our approach let us take x = x< € (0,1). Then T(h ) satisfies (3.5). Next, if 
0< 2,_1 >(O), ^(2*-i)(i) t an( j [<£(2*-i)(x+) - ^'-^(x-)], s = 1, are available, then the 

numerical quadrature formula 

T r (h) = T{h) - 

+ E (§yr h 1 * (3.6) 

will satisfy /„* <f>(t)dt - T p (h) — E p {h\x ) = 0(h 2p ) as h ->■ 0. Obviously, T p (h) with p > 2 can not 

be used as part of a quadrature method for integral equations as in (1.8), since & F(t , y{t)) and 
hence <£ (fc) (t), k > 1, are not known. We, therefore, modify T p (h) by using suitable approximations 
for the <f>W(t). Below we illustrate this approach in detail for p = 2. 

3.2.1 Modification of T 2 (h) 

We again start by taking x = x, € (0, 1). Letting p = 2 in (3.6), we thus have 

T,{h) = T(h) - ^ {(f (1) - «0)] - [«x+ ) - 4>'i.x-)}) h\ (3.7) 

and, therefore, J ^ <l>(t)dt — T 2 (h) 0(/i 4 ) as h — >■ 0. We can maintain an error of the order of 

h 4 by approximating the quantity inside the curly brackets on the right-hand side of (3.7) with 
an error of h 2 . As we want to be able to preserve the form of the equations in (1.8), we need to 
express the relevant approximations solely in terms of the F(xj, y(xj)), j = 0, 1, ..., N. Although 
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this can be achieved in various ways, we suggest the following route that seems to be the simplest 
mathematically and also very effective computationally. 

We start by breaking up <f>'{t) in the form 

4>'(t) = 0o,i (*, t)F(t, y(t)) + g(x, t) jF{t , y(t)). (3.8) 

We compute g{x,t) and 0 M (x,t) exactly since g{x, t) is given. Thus, only £ t F{t,y{t)) remains to 
be approximated. 

Approximations to f t F(t,y{t)) at t = 0 and t = 1 are provided by the one-sided three-point 
differentiation formulas 

Q'(0) = -L[-3Q(0) + 4 Q(h) - Q(2h)] + \Q'"m 2 , 0 < £ < 2h, (3.9) 

111 O 

and 

Q'{ 1) = -t[3Q(l) - 4Q(1 - h) + <3(1 - 2 h)] + 1 - 2h < £ < 1, (3.10) 

and we use these in the approximations for <t> (0) and ^ ( 1 ) • For a detailed discussion of differenti- 
ation formulas see, e.g., Hildebrand [H]. 

As for the term [<f>'{x+) - <£'(*-)], we have 

4>'{x+) - <t>'{x~) = [00,1 (x, x+) - go,i{x,x—)]F(x, y{x)) = <5i(x)F(x, y{x)). (3.11) 

Note that £F(f, y(0)|t=x is absent from (3.11) since g{x,t) is continuous at t = x. 

Combining all the above, we obtain the “numerical quadrature formula” T 2 (h) given by 

T 2 {h) = T(h) - ^{g{x i ,l)(SF N -4F N - l + F N - 2 )-g{x i ,0){-SF 0 + 4F 1 - F 2 )} 

~ ^{ffo.i(*i,l)^Ar-l/o,i(®,-,O)Fo} + ^i(x,0^ for x - x t € (0,1), (3.12) 

where Fj = F(xj,y(xj)) for short, and we have used the fact that B 2 = 1/6. This completes the 
treatment for x = Xi £ (0, 1). 

Remark. One might think that the break-down of <t>'{t ) as in (3.8) in order to apply the differ- 
entiation formulas of (3.9) and (3.10) to ±F{t,y{t)) is redundant, and that these formulas can be 
directly applied to While this is true for x = x,-, i = 2, 3, ..., N- 2, it fails to be true for x = x x 

and x = Xiv-i- The reason for this failure is that when x = x x = h or x = x N -i = 1 - h, g 0 , i(x, x) 
does not exist, hence <f>(t) is not differentiable on (0,2/i) or (1 - 2h, 1), respectively. Thus, the 
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approximations to <£'(0) and <j>'{ 1) by (3.9) and (3.10), respectively, cannot have errors of the order 
of h 2 . (Actually, the errors are 0(1) as h 0, at best.) Finally, the simplicity of the correction 
term in (3.11) coming from the point t = x is also a consequence of (3.8). 


When x = x 0 = 0 and x = Xjy = 1, the integrand is M times continuously differentiable 
for tel, hence 4>'(x+) - 4>'{x - ) = 0 in (3.7). Consequently, (3.12) is now modified to read 

f 2 {h) = T{h) - A{ 5 (o,i)(3F iV -4F w _ 1 +F w _ 2 )-<7(0,0)(— 3F 0 + 4F X - F 2 )} 

h 2 

- — {5o,i(0, 1)F n - 5o,i( 0,0+)F o } for x = x 0 — 0 (3.13) 

and 


T 2 (h) = T(h ) 


{<7(1, 1)(3 F n - AF n _ x + F n _ 2 ) - <7(l,0)(-3F o + 4F X - F 2 )} 
h 2 

{?o,i(l, l~)Fjv - 50 , 1 ( 1 , 0)F 0 } for x = x N = 1. (3.14) 


The “numerical quadrature formula” that is defined through (3.12)-(3.14) thus satisfies <p(t)dt— 
T 2 (h) = 0(h 4 ) as h — > 0, uniformly in the x,- and N (if M > 4). 


3.3 Modification of T p (h), p > 3 

Again let us begin by taking x = x, € (0, 1), and consider T p (h) in (3.6). It is sufficient to replace 
the coefficients of h 2> in the two summations there by approximations whose errors are of order 
h 2p ~ 2s , s = 1, 1. Then, the resulting modified T p (h), which we call T p (h), will maintain an 

error of order h 2p . We do this as follows: First, we break up 4> ( ~ 2, ~ x \t) in the form 

(3.15) 

Next, we approximate -£^F(t, y(t)), p = l , ..., 2s- 1, at t = 0 and t = 1, by one-sided (2p- 2s + ^)- 
point differentiation formulas, involving Xj, 0 < j < 2p — 2s + p — 1, when t = 0, and x,-, 
N - (2p - 2s + p) + 1 < j < N, when t = 1. All of the go, 2 s-i-n{x,t) at t = 0 and t = 1 are 
computed exactly. 

As for the term [^ 2l-1 ^(x-|-) — <^ 2j-1 )(x— )], we have from (3.15) and from the assumption that 
g(x, t) is continuous for t E /, that 

4> {2, ~ 1) {x+) - <£ (2,-1) (a:— ) = £ ^ ~ ^ 2,-1 - P (x)^F(t,y(t))\ t=x . (3.16) 
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We approximate £rF(t,y{t))\ t=x , fi = 1, ...,2s- 2, by (2p-2s + p) -point differentiation formulas 
in which x is in the center of the point set or as close to the center as possible. All of the i5 2j _i_^(x) 
are computed exactly. 

Note that all of the differentiation formulas above will have errors of order h 2p_2j , s= 
under the assumption that M > 2p - 1, as can easily be shown. 

The “numerical quadrature formula” T p (h) that is obtained by the approximation procedures 
above obviously satisfies / 0 X <f>{t)dt - T v (h) = 0(h 2p ) as h -+ 0, uniformly in the x, and N (if 
M > 2 p). 

In the next sections we shall refer to T p (h) as a numerical quadrature formula even though it is 
not one in the true sense of the expression. 

3.4 The New Quadrature Method from T 2 (/i) 

We close this section by giving the new quadrature method for (1.1). It is defined through the 
following system of equations 

N 

y,=r(x,) + h'^2"g(x i ,Xj)Fj 

j=° 

- 1)(3Fv " 4Fv " 1 + Fn ~ 2) ~ 0)(-3F ° + 4Fl ~ F2)} 

— l)-fiv — 9o,ii x ii 0)^o} + l = 1» 2, ..., iV — 1, 

N 

yo = r(x 0 ) + hJ2"d{ x o, x j)Fj 

j=0 

- A{ ff (0,l)(3F N -4F JV _ 1 + F J v- 2 )-</(0,0)(-3F o + 4F 1 -F 2 )} 

- ~{9o,i(0,1)F n - 9 oAO,0+)Fo} 

N 

y N = r(x N ) + h^2,"9i x N, x j)Fj 
i = o 

- Ya {g( 1. l)(3Fjv - 4Fjv_i + F n _ 2 ) - g( 1, 0)(-3F o + 4F X - F 2 )} 

- ^Ki(1,1-)^-5o, 1 (1,0)F o }. (3.17) 

Here F; = F(x,-,y,) and y is the approximation to y(x,). 

If we now write the system in (3.17) as 

Vi = $i(yo»yi,— ,Vn), i = 0 , l,...,iV, (3.18) 
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then the method of successive approximations takes the form 


y (0) _ y(°)(x f ), i = 0, 1, N, 

y,- m+1) = ^»(yo m) ,yi m) ,-iy/T ) ), * = o, ™ = 0,1,2,.. 


(3.19) 


4 Preliminary Results on r(/i) and T p (/t) 

4.1 Bound on Trapezoidal Rule Error 

The next theorem bounds the error in the trapezoidal rule in terms of the modulus of continuity 
of the integrand. 


Theorem 4.1 Assume that Q(t) is integrable on I , and denote by T^{h) the trapezoidal rule 
approximation to Jq Q(t)dt. Then 

|E q (/»)| = l^ 1 Q(t)dt - r 9 (A)| < u Q (h), (4.1) 

where u>q is the modulus of continuity of Q(t) on I. 


Proof. We have 


and 


/ Q(t)dt = £ / Qm 
Jo ,-o •'*. 


T Q (h ) = f; \ r we*-) +<3(*.-+i)i*- 

,=o z Jx ' 

Subtracting (4.3) from (4.2), we obtain 

E Q (h) = g i {£‘ +1 [Q(t) - Q( Xi )]dt + £' +1 [Q(t) - Q(* i+1 )]«&} . 

Taking absolute values on both sides of (4.4), we next obtain 

1 f r*«+i 'i 

l^)l<E 2 U. |Q(t)-Q(x i )|d« + 1 \Q(t)-Q(x i+1 )\dty 

The result in (4.1) now follows from (2.7) and from the fact that x< +1 — X { = h = 1/N. □ 
Our next result is an application of Theorem 4.1 with Q(t) = g(x,t)F(t,u(t)). 


(4.2) 

(4.3) 

(4.4) 


(4.5) 


Theorem 4.2 Assume that g € C(f2) and that F, Fo,i € C(A), and define G(x,t) = g(x, t)F(t, u(t)). 
Assume also that u(t) is such that u(t) € J for t 6 I, and G(x,t), as a function oft, is integrable 
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on I when x € 7. Denote by T G (h;x) the trapezoidal rule approximation for f^ G(x,t)dt. Then 

|F G (h;x)| = If' G{x,t)dt-T G {h;x) 

Mo 

< ||F|| u g (h) + Hill [u F (h) + UFo.xll «„(*)] , (4.6) 

uniformly in x £ I. Here ||< 7 || is the L^-norm of g(x,t ) on Q defined by ||(?|| = raax( Ii( ) e n |j(x,i)|- 
Proof. From Theorem 4.1 

|F g (/i;x)| < sup{|G(x,t) — G(x,t') \ :t,t' € 7 and |f — t'\ < h} . (4.7) 

Now 

G(x, t) - G{x, t') = F(t, «(*))[&(*> 0 “ 5(*i *')] + 9 ( x > 0[^(*» u (0) “ F ( t> > “(0)1 ( 4 - 8 ) 

and 

F(t, u(t)) - F(t\ u(t')) = [ F{t , u(t)) - F(t, *(0)] + W, u(t')) - F(t u(t'))] (4.9) 

and, finally, by the mean value theorem, 

F{t,u{t)) - F{t,u{t')) = Fo,!^^)^) - for some w 6 J. (4.10) 

The result now follows by taking absolute values in (4.8)-(4.10) and maximizing over 7, Q, and A. 
We leave the details to the reader. □ 

4.2 Bound on Error in f 2 (h) 

We now proceed to the corrected rules T G (h]Xi) with G(x,t) = g{ x it)F(t,u{t)) that are obtained 
from the trapezoidal rule r G (/i;x,) for /„* G(x,,t)dt of Theorem 4.2 exactly as described in the 
previous section. It is sufficient to examine the details of the case p = 2 as all other cases are 
treated in exactly the same way and the conclusions are the same for all values of p. 

From (3.12)-(3.14) it is clear that the correction to T G (/i;x,-), i = 0, 1, is of the form 

2 

T G (h; x.) - T G {h- x.) = h ]T [oyFfo, «(*,)) + «(* w ^))] 

j=o 

+ h 2 'u(xo)) “h N) u {%n)) “b 

There are two important points to be noted here: (i) The number of function values F(x j ,u{x j )) in 
this correction is fixed and is thus independent of N . (ii) The coefficients ctij, fiij, and 7tj are some 
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constant multiples of 0 (x„O), s(x,-, 1 ), 5 o,i(*», 0 ), So,i(z»,l), and <Ji(x<), and thus are uniformly 
bounded in i,j, and N if g 6 C(f2) and g 0 ,i € PC(Q). Thus, if F 6 C( A), we have for all i 

\T G (h; Xi ) - T a (h ; x,-)| < (C[ 2) h + Cf > h 2 )||F|| ( 4 . 12 ) 

for some positive constants Cj 2) and C 2 that depend on g but are independent of F, i, and h. 
Combining all this with Theorem 4.2, we have the following important convergence result on 

f G (h; Xi ). 

Theorem 4.3 Assume that g G C(fi) and that g Qyl G PC(Q). Assume also that F, F 0> i G C(A). 
ZeJ u G C(I) and u(t) G J for tel . TAen, /or i = 0, 1, ..., iV, we /mue 

|F^(/i; X|)| = I / (7(x,*, t)dt — T G (/i; x,*) < D 2 (h) + au; u (/i), (4-13) 

Mo 

where a > 0 is a constant and D 2 {h) is a function that goes to 0 monotonically as h — > 0, and both 
are independent o/i, N , and u(t ). 

The result of Theorem 4.3 is important since it states that T G {h\ x t ) — ► /q 1 G(x t , t)dt as /i — ^ 0 
uniformly in i = 0, 1, ..., N. (For a and D 2 {h) see Theorem 4.4 below.) 

4.3 Bound on Error in T p (h) for Arbitrary p 

So far we have analyzed the properties of T G (h; x,). We would now like to summarize the properties 
of Tp(h]Xi) for arbitrary p > 1. Note that T G (h;Xi) is simply the trapezoidal rule T G (h;Xi ) 
throughout. 

Theorem 4.4 Provided g G C(ft) and g 0ti G FC(fi), 1 < i < 2p - 3, and u G C(I), we can 
construct the corrected trapezoidal rule T G (A;x f ), which is of the form 

f p G (h ; x,) = T G {h- Xi) + £ f£ 2 A$h fc ) F(xj, ufo)). (4.14) 

j=0 \k = l / 

/fere fhe Aj p * depend on g but not on F, and can be bounded independently ofi, j, and N , and the 
number of the nonzero A'f-l is fixed and thus independent of N . As a result, under the conditions 
that u(t) € J for t € I, and F, F 0 i € C(A), (4-13) can be generalized to 

\E p (h; x,)| = \[ 1 G(Xi,t)dt-f?(h-, Xi ) < D p {h) + <ru u {h ), (4.15) 

Mo 
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where a > 0 is a constant and D p (h) is a function that goes to 0 monotonically as h — > 0, and both 
are independent of i, N , and u(t ). While o is the same for all p , D P (K) depends on p. Specifically , 

o — 11^0,1 1| and D P (h) = \\F\\ u g {h) + X C k p ^h k + ||^|| w F (/i), (4-16) 


where, for each k, Cjf 1 is an upper bound on that depends on g but is independent of 

F, i, and N. (Ifp= 1, then the A\f k and hence C ( k p) are all zero.) 

5 Existence and Uniqueness of Numerical Solution 

In order to solve the problem of existence and uniqueness of the numerical solution we need to 
enlarge the set J = [R u R 2 ] by an arbitrarily small amount r\ > 0 that will be fixed later. Hence, 
we define J' = [f?;, R' 2 ] = [Ri -ti,R 2 + g] and A' = /x J'. We also define 

||F 0ll |r = max |F 0 ,i(t, to)|. 

We start with the following lemma. 

Lemma 5.1 Assume that g(x,t ) is as in Theorem 44 so that we can define T°{h; x,) as in (4-U)- 


Define also 


Z p>i {h ; { w k }) = h X" g(xi , Xj)F{ Xj , Wj ) + X (X A %l h } F ( x i’ w i)- 

j- 0 j= 0 \* = 1 / 


Then } provided that F , Fo,i G C(A') and Uj , Vj E /, j — 0, 1, we have 

\Z Pii (h] {li*}) — Z P} i(h; {^})| < p(h)\\u — u|| } i = 0, 1, •••> N, 


M (*) = n^iir 


where 


where [[#] is as defined by (24) and C (p) are as described in Theorem 44, and 

||«-u|| = o m^K-t; J |. (5.4) 

Remark. Note that Z p<i (h; {u{x k )}) = T p G (/i;x,) when G{x,t) = g(x,t)F{t,u{t)) for an arbitrary 
function u(t). 

Proof. The difference W Pti = Z Pii (h; {u fc }) - Z Pi< (/i; {u fc }) can be expressed as 

W Pii = hY,"d{xi, u i) - F ( x D v i)] + X f Y, A ijl h J u j) “ F ( x n v i)l 

j - o i=0 \k = 1 / 
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which, upon applying the mean value theorem, becomes 


W Pt i — 9{ x ii x j)Fq,i { x j > w j)( u j u j) + ^2 i ^2 Fo,i{ x ji w j)( u j - 

J—0 j= 0 V — 1 / 

where Wj is between Uj and Vj for each j, and hence Wj G J\ j = 0 , 1 ,..., AT. Taking absolute 

values, and maximizing appropriately, we obtain from this 


I^M-I < 


hfy\g{x il x j )\+ f^Ci P) h k 

j=0 k = l 


WFo.iW Il u — u il- 


N 

Now h y^ ;/ |g(xj, Xj)| is the trapezoidal rule approximation for /q 1 |y(x*, t)\dt, and Theorem 4.1 
j=o 

applies to it. Invoking also the definition of [[<7] we thus have 

N 

j=0 

The result now follows. □ 

We now consider the nonlinear system 


Vi = r{xi) + Z Pii (h;{y k }), i = 0,l,...,N, (5.5) 

that results from applying the corrected trapezoidal rule T p (h;Xi) to (1.1). In the next theorem 
we show that there exists a unique solution for the y*, i = 0, 1, TV, for all large N , under suitable 
conditions. When p = 1, 7* (/i; x t ) = T°(h;Xi) and these conditions are almost the same as those 
given in Theorem 2.1. For p > 2 we should impose sufficient differentiability conditions on g{x,t) 
so that Tp can be defined, as mentioned earlier. 

Theorem 5.1 Assume that all the conditions of Theorem 2.1 concerning r(x), y(x, t), and F(t , w) 
hold. Assume , in addition , that g(x,t) satisfies the differentiability conditions of Theorem 4-4 so 
that T p G is defined . Assume , by extension , that F, F 0j i € C(A 0 ), where A 0 = I x [Ri — 770, + %] 

for some > 0 and thus A 0 D A. Let the sequence of successive approximations {y- m \ i = 
0, 1, ...,iV}5£ =0 be generated from (5.5) according to 

yj 0) = y (0) (*.)> * = 0,1,..., AT, 

yl m+1) = r(*j) + Zp,.'(h;{y* m) }), i = 0, 1, ..., N; m — Q, 1, ... . (5.6) 

(Here we recall that the function y^°\x) is the initial approximation in (2.6).) Then there exists a 
constant 77 E (0, r) Q ) and a positive integer No, such that the following hold: 
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(i) y| m) eJ' = [Ri - ??, Rt + v], i = 0, 1, .... N, m = 0, 1, for each N > N 0 . 

(ii) lim^oo y^ m) = y< 6 J i = 0, 1, N, exist for each N > N 0 . 

(Hi) yi, i = 0, 1, jV, is the unique solution to (5.5), for each N > N 0 - 

Remark. As will follow from the proof below, such an 77 can be picked and once this is done any 
smaller 77 will render the theorem valid. Thus, 77 can be picked arbitrarily small. 


Proof. We start by picking 77 € (0,770) such that fl>]||F 0 ,i|r < 1 , which is possible from the 
assumption that F 0 ,i € C(A 0 ) and from ( 2 . 5 ) in Theorem 2.1. Thus the choice of 77 is independent 
of N. With 77 thus fixed, we next pick a positive integer N, thus h = l/N, for which p(h) < 1, 
with p{h) as defined in (5.3). This is possible as lim A _ +0 /u(h) = fl>H|F’o 1 iH / by the assumption 
that g e C(Q) that implies that \g\ € C(fi) so that lim^o W| s |(/i) = 0 . Moreover, we have 
fi(h) < fi(h) < 1 for all h < h or all N > N. 

From Theorem 2.2 we have 


U7j,<m)(/i) < max{u; y (o)(/i),u; r (/i) + ||F|| u g {h)} = p(h), m — 0,1,... , (5.7) 


and p(h) — ► 0 monotonically as h — y 0. Define now 

e(h) = D p (h) + ap{h), 


(5.8) 


with D p (h) and a as in Theorem 4.4. Thus e(h) -> 0 monotonically as h ->• 0 as well. Therefore, 
there exists a positive integer Nq > N , 1/Nq = ho < h. for which 


({ho) < T}[ 1 - l*(h)] - ~(h) < TJ ' 


(5.9) 


Obviously, with this h 0 we have 

•« gww i 5 5 ^ < *, » - 0. 1. 2 <*•»> 

At this point it is worth recalling that Z p j(h; {tu(xfc)}) is the corrected trapezoidal rule T^{h\ i<) 
for the integral g{xi,t)F{t , w(t))dt. 

Let us set m = 0 in (5.6) and (2.1). Upon subtraction we obtain 

ytU - y U) (x< ) = Z P} ,(h-,{y i0) (x k )}) - f g{x u t)F{t,yW(t))dt, (5.11) 

Jo 
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and, by (4.15), (5.7), (5.8), and (5.10), this gives 


||e (I) || < e(h) < 77, (5.12) 

where we have defined 

e ( m) = y (">) _ y( m ) (a:,), i = 0, 1, ..., N, and \\e™\\ = max |e| m) |. (5.13) 

By (5.12) and the assumption that y^^x) € J for x € /, it is clear that y^eJ', i = 0,1, ...,7V. 
Let us next set m = 1 in (5.6) and (2.1). Upon subtraction we obtain 

y«' 2) - y (2) ( x >) = [ZpAh {y* 1 *}) - z P A h \ {y (1) (x fc )})] 

+ [z Pli (fc; {y (1) (x fc )» - £ g( Xi ,t)F(t,yW(t))dt\ . (5.14) 

Applying Lemma 5.1 to the first brackets, and (4.15) to the second ones, and also invoking (5.7), 
(5.8), and (5.12), and, following that, (5.10), we obtain 


||e (2) || < e(*) + / x(A)||e^>|| < e(h)[l + n(h)] < v . (5.15) 

/ r\ \ 

Therefore, y,- € J 1 , i = 0, 1, ..., N, too. Proceeding by induction, we can show in general that 

m — 1 

lk ( '" , ll < e(/i) < €(*) ^Em(A)]’ < (5.16) 


< 7=0 


as a result of which yj m) € </', i = 0, 1, ..., N, for all m. This proves part (i) of the theorem. 

For part (ii) we proceed in the standard way. From (5.6) we have 

y.- m+1) - y< m) = Z P Ah\ {yi m) » - Z p ,i(h-, {yi m " 1) », m = 1,2, ... . (5.17) 

Since y- m * € J' for all i and m, Lemma 5.1 applies and we have 

||y( m+1 )-y( m) ||</i(fe)||y( m )-y( m -i)||, to = 1,2,.... (5.18) 


The result now follows by the fact that n(h) < 1 since h < h 0 , which implies that {y- m) , i = 
0, 1, ...,iV}“_. 0 is a contractive sequence arid thus has a limit. 

The proof of part (iii) follows from the continuity of the function F on A'. □ 

It is worth pointing out to the similarity between (5.18) for the numerical solution and 


max y {m+l) {x)-y (m \x) 


< Li max 
~ xei 




(x)-y< m - 1 >(*)|, to = 1, 2 


for the analytical solution, with fi as in (2.5). 


(5.19) 
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Since we can pick 7 ? arbitrarily close to 0 and thus ||F 0 ,i|r arbitrarily close to ||F 0> i||, we see 
that for very large N, hence very small h, y,(h) in (5.3) is arbitrarily close to y. in (2.5). That is to 
say, the discrete successive approximation procedure converges practically at the same rate as the 
continuous one does. 


6 Accuracy of Numerical Solution 

With the existence and uniqueness questions resolved, we now turn to that of the accuracy of the 
numerical solution y ; , i = 0, 1, ...,N of (5.5). Our proof proceeds along the same lines as that of 
Theorem 4.2.1 in [K, pp. 114-115]. 


Theorem 6.1 Under the conditions of Theorem 5.1, the numerical solution with N > N 0 satisfies 

He|| = max |y, - „(»,) | < |E?(M<)I, (6-1) 


where G(x, t) = g(x,t)F(t, y{t)). 


Proof. Subtracting (1.1) with x = x,- from (5.5), we can write 

y,- - y(x t ) = [Z p ,i{h] { y k }) - Z p<i (h; {y(x fc )»] 

+ Z Pii {h;{y{x k )}) - f y(x,-, t)F(t, y(t))dt . (6.2) 

JO 

Since y(x) € J for x G I and y, e J' , i = 0,1, ..., N, Lemma 5.1 applies to the expression in the first 
brackets. The expression in the second brackets is nothing but E p {h',Xi), the error in T p {h,\ x,) . 
Thus, taking absolute values, we have 

M = \y, - y( Xi ) | < y(h)\\e\\ + |F p G (/i;x,)|. (6-3) 

The result in (6.1) follows by maximizing both sides of (6.3) and by using the fact that y.(h) < 1. 
We leave the details to the reader. □ 

Since for all N > N 0 we have fi(h) < n(h 0 ) < 1 and thus 1/(1 < 1/(1 -fi{h 0 )), we realize 

from Theorem 6.1 that the accuracy of the numerical approximations y* is determined strictly by 
that of the numerical quadrature formula underlying the quadrature method. 

In subsection 3.1 on Euler-Maclaurin expansions we proved that E p (h;Xi) = 0(h 2p ) as h ->■ 0, 
uniformly in i and N, under the assumption that M > 2p (p > 1). This also produces the result that 
E p (h;Xi) — 0(h 2p ) as h — » 0, uniformly in i and N, whenever M > 2p (p > 1), as we have already 
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shown. This result covers all cases except some in which M is an odd integer. In case M = 2p — 1 
(P > l)i we have E p (h ; x*) = 0{h 2p ~ l ) as h -> 0, again uniformly in t and N. (This time E p (h; x*) 
involves £)> hut is given by an integral representation involving periodic 

Bernoullian functions that we will skip.) This produces the result that E^(h;Xi) = 0(h 2p ~ l ) as 
h — > 0, uniformly in i and iV, whenever M = 2p — 1 (p > 1). This then covers all possible cases 
that can occur. Both bounds on E p (h; x,) will be of use below. 

Also note that for arbitrary M the rules T G , q = 1, ..., [(Af + 3)/2j, are all well defined. Using 
the facts mentioned in the previous paragraph, it can be shown that T g m+1 )/ 2 j and T^ Af+3 )/ 2 j 
have errors of the same order, namely, 0(h M ) as h — > 0. Thus, there is no advantage to the rule 
T G M+3 )/ 2 j when we know that y E C M (I), and we shall not consider it in what follows. 

In light of the contents of the previous two paragraphs we now discuss the various possible cases 
to which Theorem 6.1 applies. 

1. The case M = 0. Here r E (?(/), g E C(f2), and F, Fo t i £ C(A'), and no other differentiability 
properties for r, g, and F are given. From Theorem 2.1 y E C(I) only. Thus the quadrature 
rule that can be used for this case is only Xi(/i) = T(h)> namely, the trapezoidal rule itself. 
Applying Theorem 4.2, we obtain 

INI < B[ 0) u g (h) + Bi 0) u F {h) 4- B^Uyih) = o(l) as h -> 0. (6.4) 

2. The case M = 1. Here r E C^/), g E C(f2) and 01 , 0100,1 € FC(f2), and F, F 1)0 ,F 0> i E C(A'). 
From Theorem 2.1 y E C x (/) only. The quadrature rule that can be used for this case is 
Tx(h) = T(/i), the trapezoidal rule. (As we mentioned above, we disregard T 2 (h) even though 
it is well defined.) Now the error in the trapezoidal rule is of order h uniformly in i . Hence 
we have for this case 

||e|| < B[ 1] h for Tiih). (6.5) 

3. The case M — 2. Here r £ C 2 (I ), g G C(fl) and gj k G PC(Q), j + k < 2, and F and 

Fj.k, j + k < 2, are all in C(A'). From Theorem 2.1 y G C 2 (I ) only. In this case too we can 

use Ti(h) = T{h). (Again, we disregard T 2 (h) even though it is well defined.) The error in 
the trapezoidal rule now is of order h 2 uniformly in i. Thus, we have for this case 

M<B? ] h 2 for Ti(A). (6.6) 

4. The case M > 3. Here r G C M (I), g G C(fi) and g jik G PC(Q), j + k < M, and F and 

Fj,k, j + k < M, are all in C(A'). From Theorem 2.1 y G C M (I) only. In this case we can 
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use the rules T p (h), p = 1,2, [{M+ 1)/2J. (We disregard f l(M+3) / 2j even though it is well 
defined.) We then have 


ik 


B m h 2 Pj p = i } [{m — i)/2J, 
B^h M , p=[{M + 1)/2J. 


(6.7) 


Thus the maximum accuracy that can be achieved is determined by the differentiability properties 
of the exact solution y(x), which, in turn, are determined by those of r, g , and F . 


7 Applications to Two-Point Boundary Value Problems 

Example 1. Consider the two- point BVP with 

y" = 2 y 3 , 0 < * < 1, 

y(0) = 2 and j/( 1) = jj. 

The exact solution to this problem is 


General problems with ODE’s of the form y" = ay n with a > 0 and n > 1 occur in nth order 
reaction kinetics, see Aries [A]. We note that Dirichlet BC s are not the standard BC s associated 
with reaction kinetics problems (normally, y'(0) = 0 is assigned at x = 0). We use Dirichlet BC’s 
in our example, as this enables us to determine the exact solution by which we can demonstrate 
the accuracy of the corrected trapezoidal rule quadrature methods rigorously. 

We observe that, for x G [0,1] and y G [0,2], f(x,y) = 2y 3 satisfies 0 < ^ < 24 and 0 < 
f(x,y) < 8 y. Therefore, Theorem 1 in [P2] applies, and we conclude that (i) a unique solution 
y[x) € [0,2] exists, and (ii) with fc 2 = 12 and y^°\x) = r(i) in (1.4)— (1.7), y(x) = lim m _,oo y ^ ^(x) 
uniformly in [0, 1]. In turn, all the conditions of our Theorem 2.1 and hence of Theorem 5.1 as well 
are satisfied, and the quadrature method solutions via the trapezoidal rule and its modifications 
exist and are unique for all large N. 

From (5.18) and (5.3), it is clear that the contraction parameter p(h) of the sequence of suc- 
cessive approximations {2/,- m \ 1 = 0, 1, ..., iV}“_ 0 is practically the same as p, the contraction 
parameter of the sequence {j/ (m) (^)}m=o> that is S iven in ( 2 - 5 )- Consequently, we can conclude that 
the sequences {2/- m \ * = 0, 1, ..., N}m=i will converge to a prescribed accuracy in the same number 
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of iterations independently of N. We have verified this conclusion by solving the problem above 
with different values of N. 

We have applied the quadrature method via the trapezoidal rule T(h ) and also via the corrected 
trapezoidal rule T 2 (h) with N = 100 and N = 200. The results of the computations are given in 
Tables la and lb, respectively. The last columns of these tables demonstrate very clearly the orders 
of accuracy of ft 2 and ft 4 . It seems from these tables that if the quadrature formula T 3 (ft), whose 
order of accuracy is ft 6 , were used, then we would be able to achieve errors of order 10~ 12 with 
N = 100. We also note again that the computational cost per iteration of all three quadrature 
methods is practically the same, and this makes the high-accuracy methods practical. 


Example 2. In a problem concerning the analysis of heat and mass transfer in a porous catalyst, 
see Kubecek and Hlavacek [KH], the following two-point BVP is obtained: 


y" = oty exp 


( 7ft(l -y) \ 

u + /3(i -y)J’ 


0 < X < 1, 


y'(0) = 0 and y(l) = 1. 


The quantities 7,/?, and a are positive constants representative of dimensionless energy of activa- 
tion, heat evolution, and Thiele’s modulus, respectively. 

The solution y(x) can again be shown to satisfy a Fredholm integral equation of the form (1.4) 
with 




1 

k cosh k 


cosh kx sinh k( 1 — t), 0 < x < t 
sinh k( 1 — x) cosh&t, t < x < 1, 


and 

, x cosh kx 

r(l) = To ihF' 

The existence and uniqueness of a solution y(x) satisfying 0 < y(x) < 1 when 0 < x < 1 was 
proved in [P3]. There it is shown that f(x,y) = ay exp satisfies 0 < < 2ae lf3 and 

0 < f{x } y) < ae ll3 y for 0 < y < 1 , provided 7 (3 < 1 . Consequently, Theorem 2B in [P3] applies, 
and the sequence of successive approximations {y^ m ^(^)}^ =0 y^(x) = r(x) converges to the 

unique solution when we pick k 2 = ae 1 ^ . Again, all the conditions of our Theorem 2.1, and hence 
of Theorem 5.1 as well, are satisfied, and the quadrature method solutions via the trapezoidal rule 
and its modifications exist and are unique for all large N . 

We have applied the quadrature method via the corrected trapezoidal rule T 2 (ft) with N = 
50, 100, and 200. This enables us to verify numerically the ft 4 behavior of the quadrature method 
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error even though we do not have the exact solution against which to compare our approximate 
solution. If we let y^,- stand for the approximate solution at x = i/N, then, for each * = 0,1, ..., N, 
the ratios {y N<i - y?N v) / {V in ~ Van, a) should approach 2 4 = 16 as W becomes large. This was 
seen to be the case in the numerous computations that we performed. 

We have also observed that, for 7 /3 fixed, as a becomes large, the ODE has the characteristic of 
a singularly perturbed problem ey" = /(x,y), where e = 1/a and /(x,y) is independent of 6. Our 
computations suggest the existence of a boundary layer near x = 1 for large values of a. Even in 
such cases our quadrature method seems to be producing very smooth and accurate approximations 
to y(x) everywhere in [0, 1]. This is true, in particular, near x = 1 where y(x) has a boundary layer 
and hence is varying rapidly. 

In Tables 2a and 2b we give some numerical results obtained for the cases (a) a = 10 , f 3 = 
0.5, 7 _ 2 and (b) a = 100, /? = 0.5, 7 = 2, respectively. These tables show the numerical 
solution with N = 200, the differences d- 50) = |y 5 o,» - yioo,2i| and 4 100) = |yioo,2i - y2oo,4»| and their 
ratios d^ 0 ) / d\ 100) . In the absence of knowledge of the exact solution, and making the reasonable 
assumption that y2jv,2» is a better approximation than y/v,j, we can say that d\ and d i are 
almost identical to the absolute errors in yso,: and yioo,2«) respectively. In addition, since the 
corrected trapezoidal rule To (h) has error of order h 4 , d\ ^ / d\ ^ must approach 2 16 as .V 

becomes large. This is observed in the last columns of our tables. 

8 Summary and Concluding Remarks 

In this work we have considered the quadrature method solution (via the trapezoidal rule) of Fred- 
holm integral equations of the second kind described by (1.1) and (i)-(iii) in the first paragraph 
of Section 1. Exploiting the known singularity structure of the kernel function y(x,t), we have 
designed a class of corrected quadrature formulas for the integral f 0 g(x,t)F(t, u(t))dt with correc- 
tions derived from an analysis of the Euler-Maclaurin expansion associated with the trapezoidal 
rule. These new quadrature formulas allow us to improve the order of accuracy in the standard 
trapezoidal rule from N ~ 2 to N~ 2 p for arbitrary p > 2, where N+ 1 is the number of points in the 
discrete approximation. We have also shown that the accuracy of the quadrature method solution 
yi, i = 0, 1, is the same as that of the underlying numerical quadrature formula. 

One can also achieve an increase in accuracy by extrapolation provided an asymptotic expansion 
for the error involving negative powers of N is known. However, for one extrapolation the problem 
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will have to be solved for a given N and then again for 2N . The improvement will only be able to be 
achieved on the course grid at an expense that is almost 4 times that of the improved quadrature. 
With the improved quadrature we do not need to obtain another approximation with twice the 
number of points. 

Finally, we would like to comment on the solution of the nonlinear system of equations (5.5) 
for the discrete approximation y, , i = 0, 1, N, The reader is aware that throughout the paper 
we have emphasized the method of successive approximations as given in (5.6) for solving this 
system. Actually, successive approximations have served as an indispensable tool in the theoretical 
study of the new methods proposed here. In particular, the proof of Theorem 5.1 on existence and 
uniqueness of the y,, i = 0, 1, ...,7V, and that of Theorem 6.1 concerning the error in the y t rely 
entirely on successive approximations. 

In addition to being a theoretical tool, the method of successive approximations has something 
to offer as a practical numerical tool for actually solving the system in (5.5). First, it is an extremely 
easy method to implement on a computer. Next, as we have shown in the course of Sections 5 and 
6, the convergence of successive approximations in our problem has the nice property that the 
associated contraction parameter fi(h) given in (5.3) is practically independent of 77 , of A, and 
of which quadrature formula Tf is being used, since rj > 0 is arbitrarily close to zero and N is 
sufficiently large, and thus lim^o (lim/,_> 0 A*(^)) = fj, with /x as in (2.5). This implies that the 
number of iterations to reach convergence is nearly independent of N and of the accuracy of the 
quadrature formula. Consequently, with N fixed, the cost of the solution is practically the same for 
all accuracies. For these reasons the method of successive approximations may be a very efficient 
numerical tool for obtaining the y : when the contraction parameter /x is sufficiently smaller than 1, 
as it will require a small number of iterations to reach convergence. 

Of course, when the contraction parameter n is too close to unity, successive approximations 
converge very slowly and hence become quite expensive. In such a case Newton’s method may be 
very efficient as it has quadratic convergence and thus may produce the y t - with high accuracy at 
the cost of a small number of iterations. Newton’s method may be more efficient than successive 
approximations in such cases despite the fact that each of its iterations has a large computational 
cost. Now we need a reasonable initial approximation for Newton’s method in order to reduce the 
number of iterations and hence the cost. Again, successive approximations can be used to produce 
such an initial guess. Thus, this kind of a combination of successive approximations and Newton’s 
method may be a good way for determining the y, from (5.5). 
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Whatever the value of /x. we can also employ vector extrapolation methods, such as the mini- 
mal polynomial extrapolation (MPE) or the reduced rank extrapolation (RRE), to accelerate the 
convergence of the sequences of successive approximations from (5.6). As no Jacobian matrices 
need to be computed and no large scale linear systems need to be solved in applying MPE or RRE, 
this approach to the solution of the nonlinear equations in (5.5) via successive approximations and 
vector extrapolation methods may turn out to be more economical than that of Newton’s method, 
at least in some cases. For the subject of vector extrapolation methods we refer the reader, for 
example, to the review paper by Smith, Ford, and Sidi [SFS], and to Sidi [S], where a FORTRAN 
program that implements MPE and RRE in a numerically stable way is also given. More references 
to developments pertaining to MPE and RRE can be found in these two papers. 

Clearly, the problem of actual solution of (5.5) is of importance in itself and should be the 
subject of a separate publication. 

One last remark that we would like to make is that the approach of the present work can be 
applied to systems of nonlinear Fredholm integral equations, and hence to systems of nonlinear 
two-point BVP’s, almost with no modification. 
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X 

g( 10 °)(x) 

e( 20 °)(x) 

e( 10 °)(x)/e< 200 >(x) 

0.1 

4.34773 - 05 

1 . 08773-05 

3.99973 + 00 

0.2 

6.50813 - 05 

1 . 62773 - 05 

3.99973 + 00 

0.3 

7 . 67113-05 

1.91873 - 05 

4.00073 + 00 

0.4 

8 . 206 Z 3 - 05 

2.05273 - 05 

4.00073 + 00 

0.5 

8 . 21813-05 

2.05573 - 05 

4.00073 + 00 

0.6 

7 . 721 D -05 

1 . 93073 - 05 

4.00073 + 00 

0.7 

6.69113 - 05 

1.67373 - 05 

4.00073 + 00 

0.8 

5.09113 - 05 

1.27373 - 05 

4.00073 + 00 

0.9 

2.87673 - 05 

7 . 18973-06 

4.00073 + 00 


Table la. Results from quadrature method solution via the trapezoidal rule for Example 1. Here 
eO'O(x) stands for |t/j — y(xj)|, where x = x, = i/N for some i and y, is the approximate solution 

with h — l/N . 
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X 

e( 100 >(x) 

e< 200 >(x) 

e (l°°)( x )/ e (2 00 )( x ) 

0.1 

9.571jD - 09 

5.9857? - 10 

1.5997? + 01 

0.2 

1.0561?- 08 

6.6017?- 10 

1.5997? + 01 

0.3 

9.611Z5- 09 

6.0107?- 10 

1.5997? + 01 

0.4 

8.2917? - 09 

5.1847? - 10 

1.5997? + 01 

0.5 

6.956 D - 09 

4.3497? - 10 

1.5997? + 01 

0.6 

5.6527? - 09 

3.5347? - 10 

1.5997? + 01 

0.7 

4.3507? - 09 

2.7207? - 10 

1.5997? + 01 

0.8 

3.0037?- 09 

1.8777?- 10 

1.5997? + 01 

0.9 

1.5667? - 09 

9.7907?- 11 

1.5997? + 01 


Table lb. Results from quadrature method solution via the corrected trapezoidal rule f 2 (h) for 
Example 1. Here e<*)(x) stands for | y< - y(x.) |, where x = x { = i/iV for some t and y< is the 

approximate solution with h = 1/iV. 
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n 


d< 5 °)(i) 


S 50 \x)/S l00 \x) 

0.0 

2.9025077? - 02 

2.0087? - 10 

6.0447?- 11 

3.3227? + 00 

0.1 

3.1859527? - 02 

1.026Z? - 09 

9.5047? - 11 

1.0797? + 01 

0.2 

4.0907691? - 02 

2.274Z? - 09 

1.6207? - 10 

1.4037? + 01 

0.3 

5.7898657? - 02 

4. 1427? — 09 

2.7157?- 10 

1.5257? + 01 

0.4 

8.6037317? - 02 

6.8137? -09 

4.3367? - 10 

1.5717? + 01 

0.5 

1.3051017? - 01 

1.0277? -08 

6.4617?- 10 

1.5897? + 01 

0.6 

1.992036Z? - 01 

1.3817? -08 

8.6477? - 10 

1.5977? + 01 



1.5197? -08 

9.4757? - 10 

1.6037? + 01 

0.8 

4.590925Z? - 01 

9.8807? - 09 

6.1247?- 10 

1.6137? + 01 

0.9 

6.849658D - 01 

4.6387? - 09 

2.9737? - 10 

1.5607? + 01 


Table 2a. Results from quadrature method solution via the corrected trapezoidal rule T 2 {h) for 
Example 2 with a = 10, (3 = 0.5, and 7 = 2. Here t/ 2 oo(z) is the numerical solution with 
N = 200, d (50) (a;) = |y 5 o(z) - J/ioo(z)|, and d (100) (x) = |t/ioo(z) - !/2oo(z)|. 
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X 

y 2 oo ( x ) 

dW ( x ) 

d( 100 >(*) 

d( 5 °)(x)/d( 100) (x) 

0.0 

2.072705 D - 06 

6.755 D - 11 

3.36819 - 12 

2.00519 + 01 

0.1 

4.44089919 - 06 

7.612 D- 11 

4.622Z9 - 12 

1.64719 + 01 

0.2 

1. 695708 jD - 05 

1.954£> - 10 

1.24519 - 11 

1.56919 + 01 

0.3 

6.82218619 - 05 

5.196D - 10 

3.35519 - 11 

1.54919 + 01 

0.4 

2.75374019 - 04 

1.06119-09 

7.027Z9 — 11 

1.51019 + 01 

0.5 

1.111655D- 03 

1.34519 - 10 

2.26119- 11 

5.94819 + 00 

0.6 

4.485999D - 03 

1.58019 — 08 

9.37119 - 10 

1.68619 + 01 

0.7 

1.807543D - 02 

1.23119-07 

7.51119-09 

1.63919 + 01 

0.8 

7.238558D - 02 

6.28819 - 07 

3.85819 - 08 

1.63019 + 01 

0.9 

2.827760C - 01 

1.61419 — 06 

9.72119 - 08 

1.66019 + 01 


Table 2b. Results from quadrature method solution via the corrected trapezoidal rule T 2 {h) for 
Example 2 with a = 100, 0 = 0.5, and 7 = 2 . Here y 2 oo(x) is the numerical solution with 
N = 200, d (50) (x) = |s/so(*) - yioo(*)|, and d( 100 >(x) = |yioo(*) - lfeoo(*)l- 
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